clear all;clc;close all;
load('C:\CloudStation\ResearchNote\DiffusionSpin_NF\PumpTrans\PumpTrans_sendtoDW.mat')
Pin = cell2mat(inputPower);
Pout = cell2mat(outputPower);
T = regexp(TempLegendcell,'\d+(\.)?(\d+)?','match');
T = cell2mat(arrayfun(@(x) str2double(T{x}),1:7,'Uniformoutpu',false));

%% directly calcualte transmission
slope = [0.92,3.5075,5.55,8.0,10.55,12.625,15.825];
w=0.24;Prel = 1.48;
OD = slope / Prel;
abs_int_lam=zeros(length(Pin),length(slope));
for i=1:length(Pin)
    for j=1:length(slope)
        abs_int_lam(i,j) = 1.0 - 0.975*Prel*int_lamber(0.7*Pin(i,j)/Prel, slope(j), w)/(0.7*Pin(i,j));
    end
end

absorption = ((0.7*Pin - Pout)./(0.7*Pin));

figure();
hold on;
arrayfun(@(x) plot(Pin(:,x)/slope(x),absorption(:,x),'.'),3:length(T))
set(gca,'YScale', 'log','XScale','log');
arrayfun(@(x) plot(Pin(:,x)/slope(x),abs_int_lam(:,x),'-'),3:length(T))
xlabel('P_{in}/OD');ylabel('abs')
legend(TempLegendcell{3:length(T)})

figure();
hold on;
arrayfun(@(x) plot(Pin(:,x),absorption(:,x),'.'),3:length(T))
arrayfun(@(x) plot(Pin(:,x),abs_int_lam(:,x),'-'),3:length(T))
xlabel('P_{in}/OD');ylabel('abs')
legend(TempLegendcell{3:length(T)})

figure()
plot(T,OD,'-*'); 
xlabel('T');
ylabel('OD')

%% check CubeSolver function? Since we must solve I(x,y,z), CubeSolver is much slower than directly calculate transmission.

nGrid = 100; lengthD = 1.0; Prel = 1.48;
slope = [0.92,3.5075,5.55,8.0,10.55,12.625,15.825];
w=0.24;
beamshape = @(x,y,z) exp(-((x-0.5).*(x-0.5) + (y-0.5).*(y-0.5))/2/w^2);
tic
cs = CubeSolver('nGrid',nGrid,'lengthD',lengthD,'OD',slope(1),'power',0.7*Pin(1,2),'P0',Prel);
res = cs.lightPropagate(beamshape);
0.975*res.P_out
toc 
0.975*Prel*int_lamber(0.7*Pin(1,2)/Prel, slope(1), w)

%%

function res = int_lamber(Pin, slope, w)
    beamshape = @(x,y) exp(-((x-0.5).*(x-0.5) + (y-0.5).*(y-0.5))/2/w^2)/(2*pi*w^2*erf(1/(2*sqrt(2)*w))^2);
    I_xyz = @(x,y) lambertw(Pin.*beamshape(x,y).*exp(Pin.*beamshape(x,y)-slope));
    res= integral2(I_xyz,0,1,0,1);
end